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take weeks to run. For this reason, an alternative approach for solving time-varying boundary 
conditions called the transpiration boundary condition, has been adopted. 



Figure I CASE Lab Hardware 

The transpiration boundary condition requires no alteration of the initial discretized 
computational domain, thus saving the large amount of time required for rediscretization. The 
method alters only the normal vector and unsteady velocity boundary condition at each node on 
the surface of the body under investigation. The approach has a limitation however in that as 
surface deflections become larger, the validity of the transpiration assumption is reduced. In the 
present investigation outlined in section I, the transpiration boundary condition has been tested 
over a wide Mach number range (from subsonic through supersonic) for a variety of 
configurations. Results indicate that even for relatively large structural deformations, the 
transpiration boundary condition remains a very effective method for simulating changing 
aeroelastic boundary conditions. This makes the method particularly useful for defining local 
aeroelastic instabilities over a wide rang of operating conditions for geometries of practical 
interest. 

3. Alternative Supersonic Aerodynamic Prediction Methodology 

An enhancement to the STARS integrated analysis tool has been developed in order to 
improve the practicality of time-marched supersonic aeroelastic solutions in an operational 
environment. A significant drawback of the time-marching approach is the computation time 
required to complete the search. The use of time-marching solutions to the Euler equations on 
three-dimensional configurations of practical interest requires significant computational time in 
order to insure time accuracy. The goal of the present study has been to enhance the practicality 
of the time-marching supersonic aeroelastic analysis through a synergistic combination of 
aerodynamic modeling techniques and CFD in order to significantly reduce the computation time 
required. 

A variety of techniques were investigated to accomplish the task, and two methodologies 
have been implemented into the STARS code. Details of the methods may be found in Section 
II. In the first approach, a combination of piston theory with a modified Newtonian theory are 
implemented in a manner completely independent from the Euler Solver. In the latter approach, 
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the speed of the method is retained, with an ability to capture nonlinearities in the flow through a 
combination with the available steady Euler solution. A supersonic steady-state solution for the 
three-dimensional flowfield is obtained with the use of the finite-element Euler methodology. 
The flow variables obtained in the steady solution are saved, and treated as a mean flow 
condition which becomes the starting point for the aeroelastic simulation. The unsteady 
aerodynamic pressure is assumed to be a small perturbation about the mean flow at every surface 
node point through a local application of the isentropic wave equation. In this manner, 
nonlinearities due to 3-dimensional effects such as shock interactions are accounted for in the 
mean flow, and unsteady effects are treated as local perturbations. The unsteady generalized 
forces are obtained via integration at each time step, and the solution is marched in time. The 
process is repeated over the entire parameter space until a flutter boundary is defined. 

Due to the speed of the method, flutter boundary estimation may be approached through 
the use of a search procedure while time marching, as opposed to polynomial interpolation of 
damping values. The critical parameters obtained through the boundary estimate procedure are 
then used in the time marching aeroelastic solution using the finite element Euler methodology. 
Results indicate that estimates obtained by proceeding in this manner are sufficiently close 
enough to the boundary obtained by the Euler method (over the range of validity), that the CFD 
solution may then be used for refinement which requires significantly fewer transient 
calculations. 

Results obtained to date indicate that the combination of aerodynamic modeling with the 
CFD methodology can enhance the practicality of the use of time-marching CFD for supersonic 
aeroelastic calculations in an operational environment. For the configurations tested, it was 
shown that the time required to predict flutter by time marching CFD can be significantly 
reduced by first obtaining an accurate estimate using modeling techniques. Present experience 
indicates significant savings in identifying critical flight conditions for aircraft-type 
configurations. An added advantage of the presented approach is that the same grid may be used 
for the steady CFD solution, the perturbation model, and the time marching CFD solution. In 
addition, the method used can be easily implemented in any CFD algorithm capable of 
simulating supersonic flowfields with a relatively minimal increase in code size and complexity. 
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CHAPTER 1 


INTRODUCTION 


1.1. Background 

Ever increasing technological advancements have lead to the development of 
aircraft that operate at greater and greater velocities. Today, aircraft that are designed and 
manufactured to perform near and above sonic speeds are much more common. For 
example, it is typical for aircraft, from military fighter planes to modern airliners, to have 
cruise speeds in the transonic range. In the design and performance analysis of these 
vehicles, it is important to study their aeroelastic characteristics. Even for vehicles 
traveling at subsonic speeds, the flow can cause a body deformation in the structure, 
which in turn will affect the flow solution. The coupling between fluid and structural 
forces can induce unstable oscillations. Therefore, it is important to accurately predict a 
vehicle’s flutter boundaries, since inaccurate predictions could result in vehicle failure. 

Often in modem three-dimensional aeroelastic analysis, a flow solver is employed 
in conjunction with a structural solver. The solutions of the flow field and the surface 
oscillation are performed simultaneously. Many flow solvers used today employ a grid 
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and the surrounding flowfield. Movement of the body under investigation due to the 
surface loads must be accounted for. This requires that either the surface grid be modified 
to match the body deflection, or that the deflection be simulated by some means. One 
possibility for simulating the deflection is to modify the surface normals at their original 
locations. Thus, to the flow field, the surface appears to have an altered form. This 
concept is known as the transpiration boundary condition and is the focus of this study. 


1.2. Problem Statement 

In unsteady aeroelastic analysis, the position and orientation of the surface under 
investigation are a function of time. In order for the solution to be determined accurately, 
the surface deformation must be represented. The most direct representation would 
appear to be regenerating the computational grid. However, in a time stepping approach, 
this would require that the grid be regenerated at each time step, since the surface will 
change with every time step. Using this procedure with present computer capabilities, a 
solution for a single set of parameters using grid regeneration could be on the order of 
weeks. 

The present study uses an integrated computer code called STARS (STructural 
Analysis Routines) that is capable of performing the steady and unsteady fluid and 
structural analysis of flight vehicles that are required for determining flutter boundaries. 
The Euler based flow solver is a recent addition to the code and is capable of simulating 
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three-dimensional compressible inviscid flows. It uses finite element techniques with 
unstructured adapted meshes of tetrahedral elements. Flow solutions using this code are 
performed in a time-marching fashion. Thus, in unsteady aeroelastic applications, the 
surface deformation must be represented at each time step. This is accomplished usinu the 
transpiration boundary condition in which the original computational domain remains 
unaltered. The surface deflections are simulated by applying the deflected body normals at 
the undeflected body location, thus reducing the time required for determining flutter 
boundaries. 


1.3. Motivation for the Study 

In any field of study it is important to search for ways of improving 
solution characteristics, whether it be solution accuracy, expediency, or cost efficiency. A 
coupled fluid-structure time marching solution can be highly time consuming. For 
example, on present high-speed workstations, the calculation of a single fluid-structure 
transient on a three-dimensional aircraft configuration using the Euler equations may 
require over 100 cpu hours. Therefore, it is highly desirable to develop means of reducing 
the required computational time while maintaining solution accuracy. Advancements in 
computer technology are continually assisting in this task. However, faster computers are 
generally more expensive, so what is gained in speed may be lost in cost. By introducing 
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new concepts in lieu of or in conjunction with existing methods, vast improvements in 
solution characteristics can often be made without significantly adverse effects. 

Considering the time required to generate a single domain mesh for a flow 
solution, it seems highly impractical to rediscretize the computational domain at each time 
step. One alternative to this rediscretization was presented by Batina (1989) where he 
used a dynamic mesh algorithm that models a triangulated mesh as a spring network. 

Each edge of each triangle in the mesh is modeled by a spring whose stiffness for any edge 
i-j is inversely proportional to the length of the edge. 

In this algorithm, the grid points on the outer boundary of the mesh are fixed and a 
predictor-corrector procedure that iteratively solves the static equilibrium equations in the 
x- and y-directions at each time step is used to determine the displacements 5* and S vi at 
each interior node i. The method predicts displacements according to 

= 28 " - 5 -' 

5 y . =26". -5"- 

It then corrects the displacements using several Jacobi iterations of the static equilibrium 
equations using 



Z k A„ 
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where k m is the spring stiffness, and the summations are over all edges of the triangles that 
have node i as an endpoint. Finally, the new nodal coordinates are given by 


. n+I 


= x n +5 


n+1 


yr l =yf+5: 


This algorithm was found to produce good results when used in unsteady pitching and 
plunging studies, and also in the prediction of flutter boundaries. However, it is necessary 
to perform the displacement calculations at each node for every time step in both 
coordinate directions. Additional computations will be required if the algorithm is 
extended to three-dimensions. 

This number of calculations can be significantly reduced if the surface deflection 
can be simulated without having to alter the existing grid. One way to simulate body 
deflection is to apply a transpiration boundary condition at the surface. This is essentially 
done by rotating the body normals so that the new normals are in the same directions they 
would be in if the body had actually deflected. Thus the original body grid remains 
unaffected throughout the flutter investigation. 

Transpiration has been used effectively in simulating surface deformations in full 
potential solutions and steady and unsteady rigid-body applications in Euler equations. It 



is the purpose of this study to show the extent to which the transpiration boundary 
condition may be used in unsteady aeroelastic problems. 


1.4. Rationale 

When the transpiration boundary condition is employed, there is potential for great 
savings in time requirements. However, an expedient solution is useless if its accuracy is 
questionable. Therefore, showing the extent to which the transpiration boundary 
condition is valid will allow flutter investigations to be made much more quickly and with 
confidence in results when inside the transpiration limits. 


1.5. Flutter Analysis 

Pertaining to aerodynamics, flutter is the divergent oscillation of a surface resulting 
from a coupling between structural and fluid forces. The prediction of flutter boundaries 
has been the subject of a great deal of studies. Methods that have been used in many 
recent studies include full potential and Euler methods. Both methods have proven 
effective in many flutter investigations. 

Full potential methods offer accurate solutions in a wide range of applications with 
relatively low computational costs. As with Euler methods, there is the assumption of 
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inviscid flow. However, since these methods use an approximation of the Euler equations, 
they require further assumptions. One assumption in the development of the full potential 
equation is irrotational flow, which in most cases is a reasonable solution. These 
equations also do not allow for entropy changes across shocks. Thus, the existence of 
shocks, even in subsonic flow, will introduce inaccuracy to the solution. 

Euler methods produce a higher order, more accurate solution than the full 
potential methods. They are capable of accounting for viscous and entropy effects, 
therefore they can be used for studies of a much wider variety than the full potential 
methods, such as at high Mach numbers and with strong shocks. One drawback 
with these methods is that they are more computationally intensive. However, with ever 
increasing advancements in computer technology, these methods are becoming much more 
computationally affordable and widespread. 


1.6. Literature Review 
1.6.1. Transpiration Concept 

In 1958, M. J. Lighthill presented four alternatives for the treatment of 
displacement thickness. One of these alternatives was termed “method of equivalent 
sources”. The idea used in this method has today developed into what is known as 
transpiration. Rather than thickening an airfoil to account for the boundary layer, an 
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equivalent surface distribution of sources is used to “simulate” a thicker airfoil. This is 
done by modifying the normal velocity just outside the boundary layer to include 
additional outflow due to the boundary layer. 


1 .6.2. Transpiration in Full Potential Studies 

One early study in which transpiration was used in conjunction with a flutter 
solution was performed by Sankar, Malone, & Tassa [1981]. The study was performed 
using a full potential method. Although not explicitly stated, transpiration was used in 
“simulating” the first order bending of an oscillating rectangular wing in subsonic flow. 
This was done by applying the zero normal velocity boundary condition for the deflected 
surface at the undefiected wing position. Computations were also performed by applying 
the boundary condition at the actual surface. The authors reported making both 
computations to ensure that the differences were small, however, results were only 
presented for the transpiration solutions. 

Lift, moment, and phase results were presented for the transpiration computations 
at a Mach number of 0.24 and were compared to experiment and Kernel function 
solutions. It was concluded that the results compared reasonably well. It was noted that 
the simulated results more closely resembled the Kernel function than the experiment. 

This was attributed to neglecting viscous effects and the fact that the experiment used a 
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5% thick circular arc airfoil. Furthermore, it was concluded that the simulation accurately 
and reliably predicts unsteady subsonic potential flow. 

This transpiration boundary condition was further used in a study of a fighter wing 
in transonic flow [Malone, Sankar, & Sotomayer, 1984]. The study incorporated 
transpiration with the full potential equations to estimate the 1 st harmonic, real and 
imaginary components of unsteady surface pressures at eight different stations on an F-5 
fighter wing. Computed and experimental results were presented for three transonic Mach 
numbers (0.8, 0.9, and 0.95). It was concluded for this study that the results correlated 
reasonably well. It was further concluded that this method could be used in studying 
flutter behavior of fighter type wings at transonic speeds. 

Based on these validating results, Malone & Sankar [1985] used transpiration with 
full potential equations to study the unsteady pitching oscillation for the RAE wing-body 
model in transonic flow (Mach = 0.8). Unsteady surface pressure results were presented, 
however, no experimental data was available for comparison. 

Then, in 1986, a study was performed exclusively to compare results from a full 
potential method using the exact boundary condition and the transpiration boundary 
condition [Sankar, Ruo, & Malone], Results were presented for a NACA 64A006 airfoil 
with an oscillating trailing edge flap, a large aspect ratio wing in independent pitching and 
plunging, a plunging fighter wing, and a steady rectangular wing. 

The NACA 64A006 results were computed at a freestream Mach number of 0.875, 
The results showed that pressure distributions were within plottable accuracy of one 
another, and that integrated loads were within 10% of each other. The first harmonic out 
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of phase component of the surface pressure distribution was also determined for the flap 
case and compared to experiment. Results for both exact and transpiration approaches 
and experimental data were found to be in close agreement, the only appreciable difference 
being near the sonic line. 

The next configuration in the study was a large aspect ratio wing. Results were 
presented for three plunge velocities (corresponding to 1, 5, and 10 degrees steady anele 
of attack) at Mach = 0.77. At 1 degree, the results proved to be nearly identical. At the 
larger plunge values, the results were not as close but still very good (within 10%). 

Pitching results using the transpiration approach were also computed at Mach = 0.66 and 
compared to results from a vortex lattice method (for pitching rates of 1, 5, and 10 
degrees/second). In each case, the transpiration prediction was larger than that of the 
vortex lattice method, the worst case being over 30%. However, the vortex lattice 
method does not consider airfoil thickness effects and therefore underpredicts the airloads. 

Surface pressure distributions at four span locations were present for a Mach 0.95 
fighter wing undergoing constant plunging motion (1.5 degrees effective plunge velocity). 
Transpiration and exact approaches were compared with experiment. Both approaches 
were found to closely match experiment, except for slightly smeared shocks predicted by 
transpiration. This was suspected to be due to the use of a somewhat coarser grid in the 
transpiration case. 

The final study used transpiration to simulate steady viscous effects of a 
rectangular wing in transonic flow. The wing was analyzed at an angle of attack of 2.0 
degrees and a Mach number of 0.8. Surface pressure distributions were calculated, with 
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and without viscous corrections, and compared with experiment. Solution accuracy 
showed improvement when transpiration was used for viscous corrections. 

It was concluded from this study that for small amplitude motions, as in aeroelastic 
applications, the transpiration boundary condition provides accurate results. Furthermore, 
the authors cited a considerable savings in coding effort and memory requirements when 
using the transpiration boundary condition. 


1.6.3. Transpiration in Euler Studies 

The transpiration boundary condition was used with the unsteady Euler equations 
in a study of transonic flow past a fighter wing by Sankar, Malone, & Schuster [1987], 
The pitch oscillation of the fighter wing was accounted for by changing the boundary 
condition and leaving the original surface unmodified. In-phase and out-of-phase 
components of the surface pressure distribution were calculated at four span locations and 
compared with experiment. The fighter wing Mach number was 0.8, and the pitching 
amplitude was 0. 1 1 3 degrees at 40 Hz (zero mean angle of attack). Slightly higher 
suction levels were predicted by the Euler solver for the in-phase component near the 
leading edge. However, overall experimental and calculated results were found to be in 
very good agreement. It was concluded by the authors that the unsteady Euler solver with 
transpiration was robust enough and accurate enough to be used in aeroelastic studies. 
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Transpiration was again used in the unsteady Euler investigation of two pitching 
wings; a transonic (Mach 0.82) transport-type wing and a subsonic (Mach 0.7) rectangular 
wing [Ruo & Sankar, 1988], The transport -type wing was the Lockheed-Air Force- 
NASA-NLR (LANN) wing. Its oscillating pitch amplitude was 0.6 degrees at a frequency 
of 24 Hz. The rectangular wing’s oscillating pitch amplitude was 2 degrees at a frequency 
of 10 Hz. 

In-phase and out-of-phase components of the surface pressures were presented for 
calculated and experimental results at four span locations for each wing. For the LANN 
wing, there are significant differences in results at some locations. However, according to 
the authors, the experimental data at many locations is not considered reliable for this 
wing. When the unreliable data is ignored, results for both components show good 
agreement. For the rectangular wing, calculated and experimental results agree 
everywhere except near the wing tip for the in-phase component, which may be due to 
viscous effects. 

Midspan unsteady pressures were compared for both exact and transpiration 
methods. For the LANN wing, a steeper variation of the pressure near shock waves was 
predicted by transpiration, but away from shocks, the results were very similar. For the 
rectangular wing, almost identical results were predicted by both methods. 

For this study, it was concluded that the overall agreement of results was good. 
However, the issue of which boundary condition was more favorable for small-amplitude 
motions was inconclusive. 
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More recently, a study was performed using transpiration with an Euler method to 
simulate control surface deflections [Raj & Harris, 1993], The study included the trailing- 
edge flap deflections for a NACA 0012 airfoil and for an arrow- wing body configuration. 
Actual and simulated surface pressures were presented for the NACA 0012 airfoil with a 
10 degree flap deflection at Mach numbers 0.6 and 0.9. For both Mach numbers, the 
transpiration method produced results that were in very good agreement with actual 
results. 

The arrow-wing body configuration was analyzed using both a coarse and a fine 
grid, the fine grid having higher resolution around the wing. Simulated results were 
computed and compared with experimental data. Surface pressure distributions were 
presented at three span locations for 0, 4, and 8 degrees angle of attack with and without a 
flap deflection of 8.3 degrees. No experimental data was presented for the zero degree 
case with no flap deflection, but results for the fine and the coarse grids were in close 
agreement. Results from both grids matched well with experimental data for flap 
deflection at zero angle of attack. 

Some noticeable differences between computed and experimental results can be 
seen for the outermost station in the 4-degree angle of attack case, while there is better 
agreement at the inner stations. However, these differences are most likely due to 
limitations of Euler equations, such as neglecting viscous effects, rather than transpiration. 
Also, the finer grid appears to produce data near the leading edge that is slightly closer to 
the experimental. 
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In the 8-degree angle of attack case, results for the innermost station are very 
close. Differences are seen at the middle station, but they are not a result of transpiration. 
The differences are a result of Euler computations producing attached flow where 
experimental data suggests flow separation at the leading edge. The outermost station 
results are not as good as the innermost, but results there are comparable. 

Experimental and simulated lift and pitching-moment coefficients were also 
compared for the arrow-wing body case. Results were presented for 0, 8.3, and 17.7 
degree trailing-edge flap deflections at Mach 0.85. For each flap deflection, the simulated 
lift and pitching moment results were in very close agreement. 

From this study, the authors concluded that transpiration boundary condition was 
effective in estimating the changes in aerodynamic forces, moments and surface loading 
due to control-surface deflections, assuming that the Euler equations are capable of 
modeling the flow field. It was noted that transpiration was not suitable for simulating 
configurations that may produce geometric gaps. 


1.7. Objectives of the Present Study 

The STARS group, for which this study was performed, has recently added an 
unsteady Euler flow solver that uses transpiration to simulate surface deflections and 
deformations. The primary objective of this research is to employ this new code over a 
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wide Mach number range in the analysis of unsteady aeroelastic problems, and to 
document the extent to which the transpiration boundary condition is effective. 

It is evident that the transpiration boundary condition can be used effectively to a 
large extent in applications of relatively small displacements. It is also apparent that this 
boundary condition can be an effective tool in aeroelastic investigations. However, it is 
important to be aware of what circumstances will cause the boundary condition to 
introduce appreciable inaccuracies. Therefore, this study will present the rationale for 
when and why the transpiration boundary condition will give inaccurate results and to 
perform investigations to support the rationale. 

A secondary objective of this study is to determine how much influence grid 
resolution has on any given solution. This is a direct result of the primary objective 
because it is important to distinguish between inaccuracies due to transpiration and mere 
differences due to mesh sensitivity. 
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CHAPTER 2 


METHODOLOGY 


In this research effort, the transpiration concept was applied to various steady and 
unsteady cases. The preliminary results were obtain using a steady panel code that was 
modified to use transpiration for steady deflections. The majority of the results obtained 
from this study used an Euler based code that employs the transpiration boundary 
condition in steady and unsteady surface deflections and deformations. Computer codes 
were also developed in this study to produce actual surface deflections for use with the 
Euler code. 

2.1. Transpiration Development 

The idea of transpiration was first developed by Lighthill [1958] as a method of 
equivalent sources. Lighthill modified the normal velocity just outside the boundary layer 
of an airfoil through an equivalent surface distribution of sources to “simulate” a thicker 
airfoil. In this way, the effect of a boundary layer is present in the solution without 
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physical representation of a true boundary layer. Applied mathematically, the normal 
velocity w is 



where z is the distance from the surface and u is the x component of velocity, which takes 
the value U just outside the boundary layer. The first term is the original, unmodified 
normal velocity. The second term represents the normal velocity contribution from the 
boundary layer. Thus, the flow field “ sees” the effects of a boundary layer that is not 
physically present. 

In this study, transpiration is applied as a boundary condition such that geometric 
changes and motions are simulated through surface transpiration. The unsteady boundary 
condition states that the velocity component normal to the body surface must equal the 
velocity of the body surface, 


Vn = Vbodyn 


thus there can be no flow through the surface. With transpiration, this done by adding a 
velocity component normal to the tangential velocity so that the resultant velocity is at 
some angle to the original surface, as shown in Figure 2. 1 . 
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This velocity is then taken as the new tangential velocity. Thus the solution is performed 
on a seemingly different surface. 


2.2. Steady Panel Methods 

In general, panel methods use a distribution of singularity elements (e.g. sources 
and vortices) over a solution boundary to satisfy the solution boundary conditions. This is 
accomplished by using discrete singularity “panels” over the body surface (and possibly 
other areas, such as the wake). Combining the velocity potential of each singularity 
element with that of the ffeestream, the continuity equation is solved with the proper 
boundary conditions to give a unique solution for the velocity potential. 
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One of the conditions that is imposed in a panel method (as in other methods) is 
flow tangency; that is, there can be no flow through the body surface. Thus, in the 
solution of a steady panel problem, the boundary condition is 

V n = 0 

must be satisfied on each panel, where n is the panel normal. To apply transpiration to a 
panel method, essentially all that needs to be done is to rotate either the panel angle or the 
panel normal. 

The first investigation that was performed used transpiration with a 2D, steady 
panel code. This code, written by Arena [1993], employs the Smith-Hess panel method. 
Results were computed for an NACA 0012 airfoil with an actual and a simulated trailing 
edge flap deflection. In the actual solution, the coordinates of the original airfoil were 
changed at the flap location to create a 10 degree flap deflection. For the transpiration 
solution, the angle of each panel in the region of the flap is modified with respect to the 
free stream in the calculation of the influence coefficients. Thus, a transpiration velocity is 
not explicitly added, but rather implicitly added by an increase or decrease in panel normal. 


2.3. STARS Code 
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The computer code system used to generate the majority of the comparison data 
for this study was an extension of the original STARS (STructural Analysis RoutineS) 


computer program. STARS is an integrated Euler based FORTRAN code for the 
multidisciplinary analysis of flight vehicles (STARS users manual, 1995). Features of the 
system include structural analysis, computational fluid dynamics, heat transfer, and 
aeroservoelastic modules. The most recent version was written under the direction of K. 
K. Gupta (1995) at the NASA Dryden Flight Research Center primarily to support NASA 
flight operation, and research and development projects. Until recently, STARS used 
linearized aerodynamic theory for the prediction of the unsteady flowfields interacting with 
the elastic motion of an aircraft. This technique produces adequate results for a wide 
range of problems but has some significant limitations, such as only producing valid results 
for small perturbation flows. The technique also assumes simple harmonic motion which 
hinders studies of arbitrary or transient motions. 

Due to these limitations, an unsteady Euler CFD module was added which uses 
transpiration to simulate body surface motion. However, this addition was limited to 
simulations of single-degree-of-ffeedom, rigid body, simple harmonic motions. Since a 
primary responsibility of the STARS group was to calculate aeroelastic effects in support 
of fight test operations, the program was modified by the STARS group (Gupta, 1995) at 
NASA Dryden to allow calculation of arbitrary motion and aeroelastic effects. 

The code uses finite element analysis to derive the frequencies and mode shapes of 
the structure. The modal superposition method is then used for the dynamic structural 
response analysis. The recent addition to the code is a module that enables the 
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computation of unsteady aerodynamic forces employing the finite element-based structural 
and computational fluid dynamics computations. This code was validated with theoretical 
and experimental results. 


2.3.1. Flow Solver 

The flow solver portion of STARS is an Euler based code capable of simulating 
three dimensional compressible inviscid flows. It uses finite element techniques with 
unstructured adapted meshes of tetrahedral elements. The mesh is generated using an 
advancing front technique, which has the advantages of application to arbitrary shapes, 
varying grid density in the domain, and the ability of adaptive mesh generation in 
accordance with solution trend. 

The flow solver is comprised of different modules that perform mesh generation 
and flow solution. The five main modules are as follows: 

• SURFACE: generates the two dimensional front 

• VOLUME: generates the three-dimensional tetrahedral mesh 

• SETBND: sets the boundary conditions on the domain 

• EULER_STEADY: performs the steady Euler flow solution 

• EULER JLJNS TE AD Y: performs the unsteady Euler flow solution 
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There is also mesh geometry and flow visualization through the XPLOT and ZPLOT 
modules of the STARS system. Adaptive mesh techniques are also possible through the 
module REMESH. 

To begin a flow solution study, the user must define the curve components, surface 
components, curve segments, and surface regions that are necessary to describe the 
geometry of the problem. The user must also define the background mesh and the nodal 
spacing parameters. Taking these definitions as inputs, the code automatically generates 
the surface and volume meshes using the advancing front technique. 

The user must then define the boundary condition types for the curves and surfaces 
to be used by the preprocessor. Finally, the user must create a namelist file that assigns 
flow conditions and coefficients. The code then takes all of the user input files, plus the 
files it generates in constructing the mesh and preprocessing the flow information, and 
performs the steady flow solution. The unsteady flow solution is then performed using the 
steady flow solution as the initial condition. 


2.3.2. Aeroelastic and Aeroservoelastic Solver 

The structural solver performs nonlinear, CFD-based aeroelastic and 
aeroservoelastic analysis. The structural modeling uses a finite element method, thus 
creating a unified approach using finite element analysis for both the flow and the 
structural solution. The natural frequencies (cd) and modes (<j>) are computed by solving 
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Mii + Ku = 0 


where M and K are the inertial and stiffness matrices, respectively, and u is the 
displacement vector. The steady-state Euler solution is then performed using either an 
explicit or a quasi-implicit, local time stepping solution procedure that employs a residual 
smoothing strategy. The equation of motion in the frequency domain is 


where 

M 

K,C = 
q 

fa(t) - 


and 


fi(t) = 


Mq + Cq + Kq + f (t) + fj (t) = 0 

inertial matrix (=O t MO), and similarly 
stiffness and damping matrices 
displacement vector (=O t u) 

aerodynamic (CFD) load vector (=<PjpA), where p is the Euler 
pressure, A the appropriate surface area, and O a the modal vector 
pertaining to aerodynamic grid points interpolated from relevant 
structural nodes 

impulse force vector (= 0 T fj ) 


where f is a number of modes input by the user. The state-space form of the equation is 
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X = AX + b a (t) + b r (t) 


where 


X 

1 1 

-a- _q 

i i 

A 

i i 

1 

*> 

1 

ba(t) = 

0 

b,(t) = 

0 


and a time response solution of the state- space equation in an interval At (= t n+ i - t n ) is 
obtained as 

X„ I =e“X I +A- , [e“-r|b.(t.) + b I (t.)] 

The structural deformations u and velocities u are then computed and used by the CFD 
code to change the velocity boundary conditions at the solid boundary. Then a one-step 
Euler solution using a global time-stepping scheme is performed and the process is 
repeated for the required number of steps. 


2.3.3. Code Modifications 

The original unsteady Euler code that was used in this investigation performs 
surface deflections using transpiration only. If comparisons of these transpiration results 
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were to be made against actual deflection results, it would be necessary to produce a 
method of generating a deformed surface mesh. Essentially, the problem reduces to 
assembling a code that modifies a given surface as desired and a code for calculating the 
correct normals. Once the deflected surface is obtained with the proper normals, the 
steady Euler code can be used to perform the flow solution. 

The two codes generated to complete the process for obtaining the actual solution 
were FROMOD and SOLMOD. A flow chart depicting the process is presented in Figure 
2.3. 


Mode 

Shapes 


SURFACE 

— - — 


> FROMOD 



— . ^7 

VOLUME 


i 

SETBND 


i / 

SOLMOD t 

b J 

i 

STEADY 



Figure 2.2. Method of Surface Deflection Flowchart 


The flow chart shows that the original surface mesh does not need to be modified to 
deform the surface. A discussion of each code is presented in the follow sections. 
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2.3.3. 1. FROMOD 


The three basic components required for performing surface deflection are the 
original, unmodified mesh, the mode shapes of the different modes, and scaling factors 
specifying the generalized displacements. The unmodified mesh can be obtained by 
generating the mesh on the original surface. The mode shapes are obtained from the finite 
element structural solver. The final requirement is a set generalized displacements, which 
is provided by the user. 

The original surface is altered by displacing each node point in each coordinate 
direction, x, y, and z, on the surface by an amount that is determined from the mode 
shapes and their generalized displacements. For example, 

Xnew(l) Xold(l) Ax(l) 

where Ax(i) is determined by 

j=i 

where O a (i,j) is an array of mode shapes 

fi© is an array of generalized displacements, input by the user 
and n is the number of modes 
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Essentially, FROMOD is used to generate a deformed surface mesh. It reads in 
the surface definition and the generalized forces and prompts the user for the modal 
scaling factors. Using this data, the code generates an array of nodal displacements and 
adds these displacements to the original nodes, as previously described. Then the code 
writes the new surface file. One advantage of this procedure is that the connectivity of the 
surface is unchanged, so it can simply be transferred from the old surface file to the new 
one. The final requirement is a code to calculate the normals for the deflected surface. 

This code is discussed in the following section. 


23 . 3 . 2 . SOLMOD 

A second modification code was required to complete the process for determining 
actual deflection solutions. This code, called SOLMOD, was needed to calculate the 
surface normals for the new mesh. SOLMOD reads in the file containing the normals, 
corrects the normals, and then overwrites the old normals file with the new one. 

The computational process for modifying the normals already existed in the Euler 
code because it is used in the transpiration solution. In this process, the normals at each 
node in the mesh are calculated using an area weighted average of surface triangulations 
that contain that node. Since this coding was already present, it was simply combined 
with the proper read and write statements to create SOLMOD. 
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CHAPTER 3 


RESULTS 


The main objective of this study was to determine the effectiveness of the 
transpiration boundary condition in unsteady aeroelastic applications. This was to be done 
by performing investigations of various geometries over a wide range of Mach numbers. 

By doing so, this boundary condition could be employed in the prediction of flutter 
boundaries. The rationale is to discuss the instances in which transpiration will introduce 
error into the solution, and to generate comparison data to show the effectiveness of 
transpiration in various applications. The following section contains a discussion of the 
limitations. Subsequent sections present results demonstrating the capability of 
transpiration. 

The first set of results that is presented is for a steady, trailing edge flap deflection, 
as a preliminary investigation. The next two sets of results comprise the core of the 
research, the first being a 2 by 1 flat plate and the second being the AGARD 445.6 wing. 

In these cases, solutions were performed for actual deflections and compared to simulated 
deflections using transpiration. The final cases of the study address the issue of mesh 
sensitivity in computational fluid dynamic solutions. 
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3.1. Transpiration Limitations 


The idea of transpiration is essentially to alter the boundary condition when 
performing a flow solution. In using transpiration, if the surface under investigation 
deflects, the surface normals are modified to account for the deflection. Therefore, if the 
deflection is small, transpiration can very accurately predict the solution. This is because 
the relative position of one part of the surface to another will have little effect. However, 
as the deflection increases, the accuracy of the transpiration solution will decrease. Since 
transpiration only accounts for the orientation of the normals, the effect of translation 
between two points on the surface will not be accounted for in the solution. 

As an example of this translation problem, consider a wing that has first mode 
bending such that the free end is no longer in the same plane as the fixed end. When 
transpiration is used to simulate the bend, the normals along the wing will be altered in the 
wing’s original position. Thus the effect that a point on the displaced end has on points 
elsewhere will not be completely accurate, and increasing the displacement will add to the 
inaccuracy. 

One consequence of the translation effect is error introduced by intersecting 
shocks and surfaces. Consider a shock wave originating at the nose of an aircraft whose 
wing is oscillating in first mode bending. In the transpiration solution, the shock will 
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intersect the wing in its original position. However, in the actual solution the wing is 
displaced, thus the shock intersection will not be in the same location. 

Other problems such as this could result when there is internal surface translation 
of the body under investigation. As mentioned, however, in cases of small deflection, such 
as in flutter problems, transpiration can produce very accurate results. 


3.2. Steady Flap Case 

In order to gain a better understanding of the transpiration concept and to perform 
some preliminary investigations, a panel code was modified to employ transpiration. 
Actual and simulated surface pressures were computed for a N0012 airfoil at Mach 0.6 
and zero degree angle of attack with a steady 1 0 degree flap deflection. The results were 
obtained using a FORTRAN code written by Arena that employs the Smith-Hess panel 
method. Slight modifications to the code were made to obtain the transpiration results. 
Figure 3.1 shows the deflected flap, located at 80% of the chord. 
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Figure 3.1. NACA 0012 Airfoil with Deflected Flap 


Actual and simulated surface pressures were generated using the Smith-Hess code. These 
pressures were corrected for compressibility effects using Prandtl-Glauert correction. The 
corrected results are presented in Figure 3.2. 
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Figure 3.2. NACA 0012 Surface Pressure Results from Panel Code 

It can be seen from Figure 3.2 that transpiration is very effective in simulating the flap 
deflection with only slight differences occurring near the flap-airfoil intersection. 

This study was also performed by Raj using a steady Euler method that also 
employs transpiration. In that study, transpiration results were also compared against 
actual flap deflection results. Figure 3.3 shows the actual and transpiration simulated 
results obtained by both the Smith-Hess code and by Raj. 
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One of the primary cases in this investigation involved a flat plate 2 units long and 
1 unit wide. This type of plate is representative of configurations that have been used in 
well documented studies on panel flutter. As defined by Dixon , panel flutter is a self- 
excited oscillation of the external surface skin of a flight vehicle which results from the 
dynamic instability of the aerodynamic, inertia, and elastic forces of the system. 

Information gained from studies involving panel flutter has assisted in the pursuit for flight 
vehicles with increasing speeds. This investigation of panel flutter is gives a new approach 
to this well known problem, as it allows for flutter investigations in which the surface 
under consideration is never required to actually deform. 

The plate used in this study was centered on a surface 4 units long by 3 units wide. 
This plate and the surface mesh are shown in Figure 3.4. 



Figure 3.4. Surface Mesh for Plate Case 
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This figure shows that the grid resolution on and immediately around the plate is very 
high. This was necessary to generate a smooth flow solution over the plate. The 
freestream flows over the top of the plate, with both sides of the plate being at 
atmospheric pressure. The plate edges are pinned to the bottom surface. All surfaces are 
defined as walls and all edge boundaries are defined as far field. The plate is flat, but can 
be made to deform, as described by Section 2.3.3, or deformation can be simulated using 
transpiration. The deflection can be any combination of six bending modes of the plate. 

In plate flutter, the magnitude of deflection in each mode would not necessarily be the 
same. However in this study, all modal deflections were of equal magnitude for an 
arbitrary deflection. 

The first case compares actual and transpiration results for a generalized 
displacement of 0. 1 in all modes. As an illustration, each mode shape is shown separately 
in Figure 3.5. 



Figure 3.5. Six Vibrational Modes for Plate Case 
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The superposition of modes that produced the surface used in this case is presented in 
Figure 3.6. 



Figure 3.6. Deflected Plate Surface 


The resulting deflection produces a maximum deflection that is over 7% of the width of 
the plate. This case was selected as an initial attempt to determine limitations of 
transpiration. Arbitrary deflection was used to maintain generality; 0. 1 amplitudes were 
used to give a significant deflection. Surface pressures and generalized forces were 
compared for ffeestream Mach numbers 0.3, 0.95, and 3. In all cases, the actual and 
simulated surface pressures were very similar. The generalized forces were comparable, 
but did show some appreciable differences. 

Figures 3.7 and 3.8 show actual and transpiration surface pressure contours, 
respectively, for the Mach 0.3 case. 
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Figure 3.8. Pressure Contours Using Transpiration, Mach 0.3 


Figure 3.10. Pressure Profiles, Mach 0.3 





These pressure profiles show that there is very good agreement between the two results. 
The largest differences appear to be in regions where the plate is furthest displaced from 
its original position, which would be expected. Table 3.1 shows the generalized forces for 
the six modes of each case. 


Table 3.1. Generalized Forces for Mach 0.3 


Mode 

Actual 

Transpiration 

1 

-22.14 

-41.66 

2 

-87.48 

-91.38 

3 

-54.55 

-57.12 

4 

-4.649 

-2.229 

5 

17.65 

18.56 

6 

-57.16 

-58.33 


The generalized forces are very comparable except in modes 1 and 4. Large differences in 
generalized forces will be discussed at the end of this section. 

The surface pressure contours for the Mach 0.95 case are shown in Figures 3.11 
and 3.12. 
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Figure 3.12. Pressure Contours Using Transpiration, Mach 0.95 



Again, the pressure contours in each case are very similar. The pressure profiles for this 
Mach number are presented in Figure 3.13. 
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gure 3.13. Pressure Profiles, Mach 0.95 













These pressure profiles show excellent agreement between the two data sets with only 
slight differences being seen near some pressure peaks and in some regions of sharp 
pressure changes. The transpiration solution also seems quite capable of handlinu the 
formation of the two shocks. This can be seen again in the Mach profiles at the same 
locations in Figure 3.14. 
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Figure 3.14. Mach Profiles, Mach 0.95 







Notice how there are only minor difference near the shocks, and seemingly no differences 
at the shocks. The generalized forces for this case are presented in Table 3.2. 


Table 3.2. Generalized Forces for Mach 0.95 


Mode 

Actual 

Transpiration 

1 

-283.2 

-475.1 

2 

-1457 

-1483 

3 

-1500 

-1495 

4 

-61.35 

-34.72 

5 

233.5 

236.3 

6 

-869.1 

-886.9 


These forces are also comparable, with some noticeable differences again being seen in 
modes 1 and 4. 

Figures 3. 15 and 3. 16 show the surface pressure contours for the Mach 3.0 case. 





Figure 3.16. Pressure Contours Using Transpiration, Mach 3.0 
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As in the two previous cases, the pressures appear to be very much alike. Figure 3.17 
shows the pressure profiles. 


49 



Figure 3.17. Pressure Profiles, Mach 3.0 












Most regions show good results, however, some interesting differences can be seen in the 
region near the end of the plate. These differences could be due to the abrupt change in 
slope where the plate ends, since the edges of the plate are pinned and not clamped. Table 
3.3 gives the generalized forces for this case. 


Table 3.3. Generalized Forces for Mach 3.0 


Mode 

Actual 

Transpiration 

1 

2593 

2041 

2 

534.7 

253.1 

3 

-3708 

-3826 

4 

-848.3 

-809.4 

5 

-574.8 

-373.4 

6 

2248 

2194 


Here, the differences in generalized forces are seen in modes 2 and 5. 

As an illustration of one of the criterion for convergence, the residuals for this case 


are presented in Figure 3.18. 





Figure 3.18. Residual Convergence 

This shows a well converged solution with final residual values being on the order of 10' 10 . 

It is evident from these plots that transpiration very accurately simulates the 
deflection. However, this is not readily assumed when studying the generalized forces. 
Since the surface pressures match so closely, it can be deduced that generalized forces 
may be very sensitive to small differences in surface pressures. However, differences may 
be negligible if the structural spring force is large relative to the generalized forces. Also, 
since in one case there is actual plate deformation while in the other case the deflection is 
only simulated, the internal mesh of each domain will be different. Subsequently, as will 
be shown in later sections, the solution can be strongly dependent on the mesh. Therefore, 
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not all of the differences seen can necessarily be attributed to the tr ans piration boundary 
condition. 


3.4. AGARD 445.6 Wing 

The next case that was investigated was the AGARD 445.6 wing. This wing is a 
standard aeroelastic test configuration which has been investigated experimentally in the 
Langley Transonic Dynamics tunnel. Using this wing allows for transpiration comparison 
in practical application. Also, in this case, the surface is surrounded by the flow unlike the 
plate case where the flow is only on one side, thus adding to the complexity of the 
problem. The original, undeflected wing is shown in Figure 3.19. 



Figure 3.19. AGARD Wing 
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The wing was given a generalized displacement of 2 units in first mode bending and 2 
units in first mode torsion. The deflected wing is shown in Figure 3.20. 


Figure 3.20. Deflected AGARD Wing 

The severity of the deflection is apparent from this figure. It is so significant that an 
approximate -9 degree angle of attack is created at the tip of the wing where there should 
be none. In addition, the end of the wing has translated below its original position. 
Results were obtained for Mach numbers of 0.99, 1.141, and 2.0, using this deflection. 
Even with this unrealistic deformation at these Mach numbers, transpiration produces a 
very accurate simulation. 

The surface pressure contours on the lower surface of the wing for Mach 0.99 are 
shown in Figures 3.21 and 3.22. 



54 





Figure 3.22. Pressure Contours Using Transpiration, Mach 0.99 





Both results show the same regions of strong pressure gradients and both predict the 
formation of a shock near the middle of the tip of the wins. 

Actual and simulated surface pressure profiles at three stations for the Mach 0.99 
are shown in Figure 3.23. 
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Figure 3.23. Pressure Profiles, Mach 0.99 





The two data sets show very good agreement near the root of the wing. There is also 
good agreement near the tip of the wing except very close to the trailing edge, where the 
transpiration solution predicts a larger pressure peak. The largest differences are seen at 
the middle of the wing where the data is somewhat shifted near the trailing edge. 

Surface pressure contours are presented for the Mach 1.141 case in Figures 3.24 
and 3.25, again for the lower surface. 
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Figure 3.25. Pressure Contours Using Transpiration, Mach 1.141 


Overall, the results appear to be very good from studying the pressure contours. The 
pressure profiles for this case are presented in Figure 3.26. 
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Figure 3.26. Pressure Profiles, Mach 1.141 





Results near the root of the wing are extremely good, except for a slight difference in 
pressure peaks at the trailing edge. The agreement decreases with increasing distance 
from the root. This demonstrates the limitation described in Section 3.1 of transpiration 
not accounting for actual surface translation, since the wing is increasingly displaced from 
its unmodified position moving from wing root to tip. However, despite the very large 
deflection and rotation at the tip, the results there are quite good, with the largest percent 
error being less than 15%. 

The surface pressure contours on the lower surface for the Mach 2.0 case are 
shown in Figures 3.27 and 3.28. 



Figure 3.27. Actual Pressure Contours, Mach 2.0 
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Figure 3.28. Pressure Contours Using Transpiration, Mach 2.0 


Strong pressure gradients can be seen near the leading and trailing edges in both cases. 
The pressure change from leading edge to trailing edge appears to be more uniform in the 
transpiration case than in the actual deflection. However, in both cases the pressure 
gradient in this region is small, therefore small differences are somewhat magnified. The 
small differences are more evident in the surface pressure profiles shown in Figure 3.29. 


62 




63 


Transpiration 








Notice that the pressures away from the leading and trailing edges are almost identical 
near the root and middle of the wing. In fact, there is very good overall agreement for 
these two stations, with slight differences seen near the leading edge at the middle station. 
Again, the station near the tip shows the largest differences, as expected. 

These cases clearly show that transpiration is effective at simulating even relatively 
large deflections in transonic and supersonic flows. Another point in the application of 
transpiration is in flutter studies of a severly deflected body, such as this one. Rather than 
using methods such as deforming meshes or transpiration to deflect the surface and 
perform flutter investigations, it would be simpler and possibly more accurate to perform 
the steady solution on the actually deflected surface, then perform the unsteady aeroelastic 
analysis from this solution. The deflected surface could be obtained using methods in this 
study, and the flutter solution could be determined with confidence using transpiration, 
since it is known to be accurate in simulating small disturbances. 

As previously mentioned, some discussion of the sensitivity of the mesh to the flow 
solution is necessary so that all solution differences will not be falsely attributed to 
transpiration. This discussion is presented in the following section. 

3.5. Mesh Sensitivity 
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The original plate mesh was much less refined than the one that was used to obtain 
the results presented in plate section. The original mesh is shown in Figure 3.30, 



Figure 3.30. Original Plate Mesh 

When observing the solutions from this original mesh, it was apparent that the elements 
were too small and too few to make the solution smooth. Therefore, the solution was 
probably somewhat inaccurate. A new mesh was constructed to produce a more accurate 
solution. This mesh was presented in Figure 3.4. 

Pressure profiles are presented in Figure 3.3 1 for both the new and old actually 
deflected meshes for 0.1 generalized displacements at Mach 3.0. 
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Figure 3.31. Pressure Profiles for New and Old Plate Cases 






Some rather large differences are seen between the two data sets. In many areas, the 
differences are larger than differences that were found in the various actual and 
transpiration comparisons. Since the deflections and the flow conditions are the same, the 
difference must be the mesh. 

The generalized forces that were obtained for the plate using the old mesh are 
presented in Table 3.4 along with the values from the new mesh (from Table 3.3). 

Table 3.4. Generalized Forces for Old and New Meshes 


Mode 

Old Mesh 

New Mesh 

1 

2718 

2593 

2 

1274 

534.7 

3 

-2564 

-3708 

4 

-842 

-848.3 

5 

-2872 

-574.8 

6 

713 

2248 


There are significant differences in the two sets of values. Again, in some cases, the 
differences are larger than differences between actual and transpiration results. 

Obviously, the more refined the mesh, the more accurate the solution will be. 
However, a more refined mesh contains more elements and nodal points, thus requiring 
more computational effort. The investigator must determine a balance between accuracy 
and practicality. Due to this mesh sensitivity, a final case was investigated where only one 
mesh was used to obtain comparison results. This case is presented in the following 
section. 
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Figure 3.32. Deflected Plate, 0.5 Generalized Displacements 


The remarkable amount of deflection is obvious from the figure. The deflection is so 
significant that it produces maximum displacement that is greater than 35% of the width. 
Solutions were obtained for Mach numbers 0.3, 0.8, 0.95, and 3. In each case, a well 
converged deflected solution was used as the starting point for the transpiration solution. 
In all cases, transpiration was effective in “removing” the deflection. However, results for 
the Mach 0.95 case were not as good as the others after the same number of iterations. 

A summary of the final results for this case is given in Table 3.5. 
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Table 3.5. Summary of Results for 0.5 Generalized Displacement Plate 



Mach 0.3 

Mach 0.8 

Mach 0.95 

Mach 3.0 

Mach Number 

0.3 + 0.00001 

0.8 ±0.00005 

0.95 ± 0.04 

3.0 ±0.001 

Pressure Order of 
Magnitude 

~ 10*° 

io- 4 

10' J 

10* 

Average Percent 
Difference in GFs 

2.98e-6 

1.62e-5 

8.27e-3 

5.23e-6 


From the table it can be seen that the overall Mach number becomes approximately the 
freestream Mach number in each case. Also, the surface pressures in each case are 
essentially zero with respect to their original deflected values, which were on the order of 
10' 1 to 10°. The average percent error in generalized forces was calculated by averaging 
all of the ratios of final force value and initial force value for each mode. As an illustration 
of the generalized forces approach to zero, a plot of the generalized force time history for 
Mach 3 is presented in Figure 3.33. 
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Figure 3.33. Time History of Generalized Forces 


This plot is representative of the time history of generalized force for each of the other 
Mach numbers as well. 

A plot of average percent error versus Mach number is presented in Figure 3.34. 


71 




Figure 3.34. Average Percent Error in Generalized Forces 

This plot suggests that transpiration was less effective in removing the deflection in the 
transonic range. However, since transpiration did not appear to be less effective in the 
other cases in the transonic range, a generalization cannot be made. Overall results for 
this case are very good, considering the magnitude of the deflection that is being 
“removed” by transpiration. 
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CHAPTER 4 


CONCLUSIONS AND RECOMMENDATIONS 


4.1. Conclusions 

The primary objective of this research project was to examine the effectiveness of 
transpiration for simulation of structural deformations in steady and unsteady aeroelastic 
applications. The majority of the investigations were performed using a recently modified 
version of a highly integrated, finite element-based code for the multidisciplinary analysis 
of flight vehicles. A supplement to this code was developed in this study which allows for 
the generation of deflected meshes using modal superposition. This research 
demonstrated that the transpiration boundary condition has strong potential for 
applications in unsteady aeroelastic analysis, such as in the prediction of flutter boundaries. 
The following conclusions were reached during this investigation: 

1. The transpiration boundaiy condition is effective in simulating even relatively large 
displacements over a wide range of Mach numbers. Some the results support the 
rationale presented for when transpiration will lose accuracy. However, there is no 
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strict criteria for when the transpiration boundary condition will breakdown. 

These results show that for applications similar to the ones presented, such as 
flutter prediction, transpiration can be a very effective tool in simplifying the 
analysis. 

2. Solutions involving the application of a domain mesh can be sensitive to the 
refinement of the mesh. The researcher performing studies using domain meshes 
should investigate the sensitivity of the mesh to his or her particular application, 
and depending on the desired accuracy, employ the mesh that is the most practical. 

3. The codes developed in this study, FROMOD and SOLMOD, can accurately and 
simply be used to perform surface deflections. In cases where a surface is 
significantly deflected from its original position before it begins to oscillate, it may 
be effective to use a surface deflection scheme such as this to deform the body and 
then perform the flutter investigation from this initial condition. 


4.2. Recommendations 

The cases of this study used relatively simple geometries, as compared to full body 
configurations. According to the transpiration concept, there should be little or no 
accuracy lost in using transpiration in more complex cases. However, in instances such as 
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intersecting shocks and surfaces, it would be interesting to see how the transpiration 
boundary condition performs. 

Also, considering the findings on mesh sensitivity, it is recommended that any 
future studies using the transpiration boundary condition to compare simulated deflections 
with actual deflections begin with a documentation of the sensitivity of the solution to the 
mesh. 

Finally, since the majority of these results compared computer simulations with 
one another, it may be advisable to compare with experimental results when available to 
confirm the results from the simulations. 
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NOMENCLATURE 


P OT , p«, = velocity of sound, pressure, density, and mach number respectively, in free stream 

conditions 

ao, P 0 , p 0 , = velocity of sound, pressure, density, and mach number respectively, in mean flow 

M 0 

C p = coefficient of pressure 

C L = coefficient of lift 

y =1.4, ratio of specific heats 

5 = impact angle with free-stream velocity 

6 0 = mean flow angle 

0' = perturbed angle 

M ns = mach number normal to shock 

n = outward normal 

n' = outward normal perturbed from mean flow 

P = Pressure 

P' = pressure perturbed from mean flow 

q = generalized displacement 

w(t) = piston velocity as a function of time 

V b = velocity of body 

V s = local steady velocity 
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1. INTRODUCTION 


A prediction of aircraft flight dynamics and aeroelastic characteristics such as flutter are 
crucial to the design of modern aircraft as well as to flight test operations. Recently, there has 
been resurgent interest in the development of hypersonic vehicles which operate in speeds where 
flutter may occur during at least part of the vehicle s operation. Examples of such vehicles 
include the Space Shuttle, Pegasus, the National Aerospace Plane (NASP), and the X-34. Also 
presently under development is the X-33, which is a joint effort between NASA and major 
aerospace companies, and will be the world’s first single stage to orbit vehicle. 

Using a recently developed STARS 2 capability for aeroelastic analysis, a time-marching 
approach based on the unsteady Euler equations may be utilized to predict flutter boundaries over 
a wide Mach number range for complex three-dimensional geometries. Determination of the 
flutter boundaries is presently achieved by searching over the flight regime for potential 
crossovers between stable and divergent time history oscillations based on modal damping terms. 
This analysis is followed by interpolation of these results to determine the point at which the 

system is neutrally stable. 

STARS which stands for “Structural Analysis Routines” was developed by Gupta at the 
NASA Dryden Flight Research Center. STARS is a highly integrated code for multidisciplinary 
analysis of flight vehicles including static and dynamic structural analysis, computational fluid 
dynamics, heat transfer, and aeroservoelasticity capabilities. With the recently developed 
capability for Aeroelastic analysis using a time-marching approach based on the unsteady Euler 
equations, the aforementioned prediction of flutter boundaries may be obtained for a wide variety 
of flight conditions and geometries. 
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Due to the potentially large domain required to ensure sufficient grid resolution of a given 
geometry, however, there lies a critical drawback of the time-marching approach. Dowell 
points out that the computational time required will be on the order of P * Tp where P is the 
number of parameter combinations required and Tp is the time required for a simultaneous fluid- 
structure time marching calculation to complete a transient. In order to ensure time accuracy and 
sufficient grid resolution, the use of time-marching solutions to the Euler equations on a three- 
dimensional configuration requires a significant amount of computation time. As an illustration, 
on a present day high speed workstation, utilization of the unsteady Euler equations to calculate a 
single fluid structure transient on a three-dimensional system may demand well in excess of one- 
hundred CPU hours. Hence, identification of the flutter boundaries over the full domain will 
require many times this number. 
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2. FOCUS OF THE PRESENT INVESTIGATION 

The focus of the present investigation is to enhance the practicality of the time-marching 
aeroelastic analysis by combining aerodynamic modeling techniques and Computational Fluid 
Dynamics in order to significantly reduce the computational time required for an aerodynamic 
solution of a given geometry in the supersonic flight regime. Moreover, the error of such a 
method must be well within acceptable limits as compared with experimental data or the 
predictions of the unsteady Euler equations previously discussed. 
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3. METHODS CONSIDERED 


Since it is the CFD which requires the overwhelming proportion of the computation time, 
a method which gives an accurate estimate to later be refined is researched. Different areas of 
importance in the determination of such a method are as follows. Since this procedure will be 
implemented into STARS, ease of implementation and compatibility must be considered. Also, 
this method must be accurate over a wide range of geometric shapes and flow regimes, and not 
just limited to special cases. Simplicity is also a plus, hence a basic method would augment ease 
of rigorous testing and validation over the widest range possible. Different methods with these 
goals in mind are discussed. 

3. 1. Piston Method 

The piston method 4 is a popular technique for supersonic and hypersonic flutter 
prediction. In the form of the unsteady wave equation, this method requires only the outward 
normal vector on a given surface to determine the pressure at any point. Figure 3-1 illustrates the 
unsteady wave equation is as follows: 



Figure 3-1 Piston motion in a one-dimensional channel. 
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Due to its simplicity and ease of use, it is an attractive technique for approximating the 

surface pressure in a supersonic flow. The unsteady wave equation, however, does not take into 

j 

account the losses across a shock nor does it accurately predict pressure in a shock interaction 
area. Furthermore, since this theory is based on a point function, i.e., the pressures are only 

dependent upon local conditions, it over predicts the pressure on a three-dimensional geometry 

i 

such as a cone. This is due to the three-dimensional relaxing effect for which the piston theory 
cannot account. Also, for a relatively blunt surface with respect to the flow, the piston method 
again over predicts the pressure. As an example, in figure 3-2, the unsteady wave equation is 
used to show the pressure coefficient versus mach number for a wedge and a cone at half angles 
of ten degrees. Since the piston theory does not differentiate between the two geometries and is 
i based on the assumption that the flow is a point function, only one curve represents both 

geometries. This curve is compared with data taken from literature for a cone and wedge at the 
noted half angle. 
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Mach # 

Figure 3-2 Pressure coefficient versus Mach number for a wedge and cone at half angles of ten degrees as 

applied to Piston Theory. 

Figure 3-3 shows the percent error of the application of Piston theory to both the wedge and a 
cone as compared with data from the literature. 



Figure 3-3 Percent error of Piston theory application to the wedge and cone. 
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Clearly the unsteady wave equation predicts pressure well over a predominately two dimensional 
flow surface yet over predicts pressure on a more three-dimensional flow surface such as the 

cone with the same half angle. 

Due to the limited application range and accurate prediction of pressure via the Piston 
theory, a method which will better predict surface pressure about a more three-dimensional flow 
is considered. 


3.2. Tangent-Cone Method 

The Tangent-Cone 7 method is an approximate technique for the prediction of three 
dimensional pressure effects on a cone surface analogous to the nose of a fuselage for example. 
The equations for the surface pressure are as follows: 

48- Mi -sin 2 8 
Cp “ 23- Mi-5 

where: M„ = (0.87 • Af_ - 0.544) • sin <5 + 0.53 

C p is based on the Modified Newtonian 8 theory which yields a result as a function only of the 
impact angle 5: 

C p = K ■ sin 2 5 

Here K 8 is equal to the stagnation pressure coefficient and M ns is the mach number normal to the 
shock which is empirically determined by Edwards 9 . Results and further explanation of this 
method may be found in reference 8. 
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The Tangent-Cone method under predicts pressure for predominately two-dimensional 
flow and again over predicts pressure for flows about blunt surfaces such as the leading edge of a 
wing. Furthermore, as in the Piston Method, the Tangent-Cone method neglects the effects of 
the mean flow and any shock interactions. 

3.3. Modified Newtonian Impact Method 

In contrast to the unsteady wave equation and the Tangent-Cone procedure, the Modified 
Newtonian Impact 8 method is used predominately for very blunt surfaces relative to the free- 
stream flow or at very high mach numbers in the hypersonic range. The equation for the 
Modified Newtonian Impact theory is previously discussed with the Tangent-Cone method and is 
also listed as follows: 

C p =K-sin 2 8 

For true Newtonian flow, the parameter K is given as 2, and for the Modified theory, K is taken 
as the stagnation C p . 

Since however, this method has a limited application of mainly blunt surfaces at high 
supersonic to hypersonic speeds, other methods are again considered. 



4. Other Methods 


There are a number of other methods for prediction of surface pressure in supersonic and 
hypersonic flow. The Modified Newtonian Plus Prandtl-Meyer Method for example is another 
blunt body technique based on the analysis presented by Kaufman . Or even the Van Dyke 
Unified Method 8 useful for thin profile shapes. These methods are usually applied to specific 
cases and tend to be too specific for the task at hand. 
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5. 


METHODOLOGY 


After careful examination and analysis of several different methods to augment the 
efficacy of the STARS code in flutter boundary prediction, the Piston method is chosen not only 
to be implemented into to STARS but also to be utilized as a perturbation to an already existing 
mean flow solution, hence the Piston Perturbation Method. A brief illustration of this idea, 
utilizing a cone similar to the one used for the generation of figure 3-3, is shown in figure 5-1 
where the unsteady wave equation is modified to predict the pressure for a perturbations about 
the mean flow from 10° to 12.5°. The modified unsteady wave equation is as follows: 

i— 

p' p T v — 1 y~ l 

— = -£- 1 + ^——M 0 Sin(6-6 0 ) 

P^ PA 2 J 

where figure 5-1 illustrates this equation. 

Pco 

Figure 5-1 An illustration depicting the modified unsteady wave equation for a cone. 
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Figure 5-2 Pressure coefficient of a cone perturbed about the mean flow of 10° to 12.5°. 


Previously in figure 3-2, Piston theory was applied to a cone without taking into account the 
mean flow conditions. Figure 5-3 shows the percent error of figure 5-2, which comparatively 
shows better results than the previous situation. 


8 T 


% Error Versus Mach # Perturbed From 10° to: 




12 . 5 ° 


- □ 


% 


2 3 4 Mach # 5 6 7 8 

Figure 5-3 Percent error versus Mach number for Cp of a cone perturbed from 10° to 12.5 C 


Clearly, the percent error is less for the perturbed case than the unperturbed case in figure 3-3. 
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The Perturbation method has been implemented in STARS in the following manner: 
Initially, a supersonic steady-state solution for a three-dimensional flow field is obtained with the 
use of the finite-element Euler methodology. The flow variables, which take into account the 
non-linearities due to three-dimensional effects such as shock interactions, are saved and used as 
a mean flow to the starting point of the aeroelastic simulation. Once these variables are obtained, 
they may be saved and used for future simulations given the structure is not significantly altered. 
With the steady state Euler solution given as the mean flow, modal superposition is applied in the 
simulated aerodynamic structure to represent the aeroelastic effects. The modes of vibration are 
perturbed which represents a perturbation to the mean flow. Next, an application of the 
isentropic wave equation previously discussed is locally applied as a perturbation to the mean 
flow at every point. The equations governing this method are shown here with figure 5-4: 



Figure 5-4 Simplistic illustration of a locally applied perturbation to the mean flow. 
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with: 


a, 


P/ 

7p~ 


and A u* =V X • n' + V h -n' . 

These pressures are then used in the coupled structural dynamic solutions which are numerically 
integrated to find a generalized displacement q shown as: 

[M]-q + [K]-q = [P\ 

where [M] and [K] are the mass matrix, and the stiffness matrix respectively, obtained in STARS 
by a finite element structural analysis routine given the structural properties of the system and [P] 
is the force matrix obtained by the piston perturbation method. Once the generalized 
displacement for every point on the aeroelastic surface is determined, the values are multiplied by 
the mode shapes to determine the actual displacement for calculation of a new set of transient 
data. 

Given the new structural deflections based on the previous aerodynamic pressures, a new 
outward normal velocity is given for the next calculation of the pressures. This process repeats 
until enough transient data is acquired to determine the stability characteristics of the perturbed 
aerodynamic system. Finally, the flutter boundaries are then verified and refined with the non- 
linear time marching Euler method. 
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6. RESULTS 

It will be shown that the Piston Perturbation method accurately predicts surface pressure 
of a number of different geometrical configurations in the supersonic flow region. It will also be 
shown that the estimated flutter boundaries of a Generic Hypersonic Vehicle (GHV) obtained by 
this method are sufficiently close to the flutter boundaries obtained by the Euler method, that the 
estimated CFD solution may then be used for refinement which requires significantly fewer 
transient calculations. This GHV will also be used to compare predicted flutter boundaries, and 
run time between the perturbation solution and the Euler solution. 

To show the simplicity and accuracy of the piston perturbation solution, a simple 
perturbed wedge exemplifies the validly of this method over a wide range of mach numbers for 
perturbations about the mean flow. Figure 6-1 shows these results. 



Mach # 

Figure 6-1 A simple perturbed wedge in compression. 


14 




The Perturbation method also predicts pressure for an expansive perturbation as shown in figure 
6 - 2 . 


Wedge Cp Versus Mach # Perturbed From 12.5° 



Figure 6-2 A simple perturbed wedge in expansion. 


The errors are minimal given a perturbation approximately less than 30% of the mean 
flow conditions and w/a„ between -1 and +1 where negative and positive values represent 
expansion and compression respectively. The term w/a„ is similar to w/a J° in the unsteady wave 
equation previously noted except it is perturbed about the mean flow not the free stream. Since 
the majority of flutter analyses deals mostly with small perturbations this criteria holds for a wide 
range of mach numbers and geometries. 

The P.S. also predicts pressure about more three dimensional surfaces, such as cones, 
over a wide range of mach of mach numbers in compression as well as in expansion shown in 
figures 6-3, and 6-4 as follows. 
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Figure 6-3 A simple perturbed cone in compression. 


Cone Cp Versus Mach # Perturbed From 12.5° to: 



Figure 6-4 A simple perturbed cone in expansion. 
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Again, the P.S gives accurate results, in the reasonable situations, of the surface pressure since it 
is applied as a small perturbation to the mean flow where the three-dimensional relaxation effects 
are accounted for. 

Another advantage of the Perturbation method is it s ability to capture pressure 
predictions when nonlinearities such as shock interactions exist in the mean flow. As shown in 
the following image (figure 6-5), at M«, = 2.2, a shock induced by a rigid wedge in a fixed 
rectangular duct interacts with an elastic clamped flat plate with differential pressure on the outer 
part of the duct. 



Figure 6-5 A fixed rectangular duct with an elastically flexible clamped flat plate. 


A side view of this geometry (figure 6-6) shows the pressure contours generated by steady Euler 
equations. Notice the heavy shock interactions due to the wedge and rigid boundaries of the 

geometry. 
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Figure 6-6 A side view showing the pressure contours of the heavy shock interactions with the elastic plate. 


A magnified view for the steady state deformation of the elastic plate generated by the unsteady 
Euler equations is shown in figure 6-7. 



Figure 6-7 Steady state deformation of the elastic plate generated by the unsteady Euler equations. 


This deformation was generated by first running a steady Euler solution to obtain the 
aerodynamic properties throughout the duct. Next, the process is restarted, this time allowing the 
plate to deform due to the resulting aerodynamic forces just acquired. After the plate’s 
oscillation damps out, steady state deflection is achieved. Notice the plate s outward 

deformation due to the shock induced by the rigid wedge. 

As a comparison, pressure using the piston theory is calculated and shown below in figure 

6 - 8 . 
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Figure 6-8 Steady state deformation of the elastic plate generated by piston theory. 


Since this theory only takes into account the local conditions, the pressure induced by the shock 
is totally neglected causing large, error in the prediction of the steady state deformation. The 
Perturbation solution is now applied with the following results in figure 6-9, 



Figure 6-9 Steady state deformation of the elastic plate generated by the Perturbation method. 


Since the perturbation method perturbs about the mean flow, the pressure induced by the shock is 
accounted for giving more accurate results as compared with the Euler steady state deformation 

in figure 6-7 above. 

As another example of the capability of the Perturbation method, the pressure at three 
sectional cuts on a Generic Hypersonic Vehicle (GHV) at Mach 2.2 is calculated at a 5° angle of 
attack using the P.S. The piston method, and a steady Euler solution is also calculated for 
comparison. For the perturbation method, the GHV is perturbed 1“ about a steady Euler solution 
at 4° with the following results shown in figure 6-10. 
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Figure 6-10 Pressure comparisons at 5° angle of attack using Euler, Piston, and Perturbation method at three 

sectional cuts on a GHV. 
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As another illustration of the perturbation method a surface mesh of the baseline 


configuration for the GHV is generated (fig. 6-11). 



Figure 6-11 GHV Baseline surface Mesh 


Flutter analysis is then performed in the following manner: Given the surface mesh, a finite- 
element structural model is then developed to obtain the structural mode shapes and frequencies. 
Details of this structural modeling may be found in references 2,11, and 12. Next, a steady 
solution to the flow field at Mach 2.2 was obtained using the finite-element Euler methodology 
which is used as the mean flow condition about which the perturbation solution is applied. The 
aeroelastic simulation consists of a 9 mode solution for 705 time steps (which is approximately 7 
cycles of mode 1). The Euler solver is run over a range of 4 dynamic pressures, and the flutter 
boundary is estimated through polynomial interpolation. For the purpose of run-time 
comparison, the same procedure is used for the perturbation solution. 
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As seen in Figure 6-12, the difference among flutter boundary estimates between the two 
codes for this case is found to be less than 4%, however the difference in run times is extremely 
significant. 




CPU Time (Minutes) 


Figure 6-12 Flutter boundary and run-time comparisons. 

M=2.2, 705 time steps / transient 

Run time for the perturbation solution (P.S.) estimate is 3 minutes for each transient or 
approximately 15 minutes to identify the flutter boundary. On the other hand, run time for the 
Euler solution used to define the flutter boundary is 117 hours for each transient and 
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approximately 469 hours to identify the boundary. All simulations are run on an IBM RS6000 
3BT workstation. 

As a final illustration of the perturbation method, the following geometry is considered to 
examine the flow characteristics about a three dimensional flow field in which a highly three 
dimensional flowfield interacts with a swept wing.. The analysis is conducted in a manner 
similar to that used for the GHV case. The baseline geometry may be seen in figure 6-13. 



Figure 6-13 Test case baseline mesh. 

A wing with a biconvex cross section is swept 42 degrees which makes the design Mach number 
when placed on the 15° cone approximately 1.6. Since steep gradients require a fine mesh which 
requires more computation time, the cone is tapered aft of the swept wing. Since the taper is 
gradual and the flow is supersonic (flow effects may not propagate forward), this area may be 
coarsely meshed aiding in computational efficiency. Mach 1.6 is one Mach number considered 
along with Mach 1.3 and Mach 2.4. Figure 6-14 shows the pressure contours at mach 1.6 of this 

configuration. 
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Figure 6-14 Pressure contours at Mach 1.6 


The perturbation method and the Piston method are used to determine the flutter boundary at 
Mach 1.6. Next, these boundaries are compared with the flutter boundary determined by the 
unsteady Euler equations. Comparisons of analyses show a less than 5% error in the prediction 
of the flutter boundary between the perturbation solution and the Euler solution. The piston-only 
solution is in error by 37.5% 

The same analysis is performed at mach 2.4. Note how the shock produced by the cone 
interacts with the wing at about mid-span (fig. 6-15). This is a typical example where the 
Perturbation method is able to predict three-dimensional flow characteristics by accounting for 
the nonlinearities and three-dimensional effects in the mean flow. A plot of the steady-state 
pressure contours used as the mean flow conditions for the P.S. is shown in figure 6-15. 
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Figure 6-15 Pressure contours at Mach 2.4. 


Flutter prediction results are again similar to those seen at other Mach numbers. Differences 
between the perturbation solution and the Euler solver are again less than 5%, while the piston- 
only solution is in error by 20%. 

7. CONCLUSIONS 

Obtained results reveal that the combination of aerodynamic modeling with CFD 
methodology enhance the practicality of the use of time-marching CFD for supersonic aeroelastic 
calculations in an operational environment. This method was tested over a wide range of mach 
numbers and proven valid in all reasonable cases. Furthermore, the application of this method in 
areas of shock interactions was studied and also proven accurate. It was shown for the 
configurations tested that the time required to predict flutter by time marching CFD can 
significantly be reduced by first obtaining an accurate estimate using modeling techniques. 
Furthermore, present research also indicates significant time savings in the identification of 
critical flight conditions for aircraft-type configurations in addition to the ones discussed here. It 



was demonstrated that for highly-three dimensional configurations, and for configurations for 
which significant nonlinearities exist in the mean flow, the perturbation method developed, is 
superior to methods which cannot account for the mean flow. 

An added advantage of the presented approach is that the same grid may be used for the 
steady CFD solution, the perturbation model, and the time marching CFD solution. In addition, 
the method used can be easily implemented in any CFD algorithm capable of simulating 
supersonic flow fields with relatively minimal increase in code size and complexity . 
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